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ABSTRACT 

We present new high resolution numerical simulations of the interstellar medium (ISM) in a central 
. i? < 32 parsecs region around a supermassive black hole (1.3 x l{f Mq) at a galactic center. Three- 

' dimensional hydrodynamic modeling of the ISM (Wada & Norman 2002) with the nuclear starburst 

now includes tracking of the formation of molecular hydrogen (H2 ) out of the neutral hydrogen phase 
as a function of the evolving ambient ISM conditions with a finer spatial resolution (0.125 pc). In 
a quasi equilibrium state, mass fraction of II2 is about 0.4 (total II2 mass is ~ 1.5 x IQ^Mq) of the 
total gas mass for the uniform far UV (FUV) with Gq = 10 in Habing unit. As shown in the previous 
model, the gas forms an inhomogeneous disk, whose scale-height becomes larger in the outer region. 
II2 forms a thin nuclear disk in the inner ~ 5 pc, which is surrounded by molecular clouds swelled 
up toward h < IQ pc. The velocity field of the disk is highly turbulent in the torus region, whose 
Q . velocity dispersion is ~ 20 km s^^ on average. Average supernova rate (SNR) of ~ 5 x 10~^yr~-'^ 

' is large enough to energize these structures. Gas column densities toward the nucleus larger than 

C/3 \ 10^^ cm~^ are observed if the viewing angle is smaller than 9^ ~ 50° from the edge-on. However, 

i the column densities are distributed over almost two orders of magnitude around the average for any 

given viewing angle due to the clumpy nature of the torus. For a stronger FUV (Gq = 100), the total 
II2 mass in an equillibrium is only slightly smaller {~ 0.35), a testimony to the strong self-shielding 
^ I nature of H2, and the molecular gas is somewhat more concentrated in a mid-plane. Other properties 

' of the ISM are not very sensitive either to the FUV intensity and the supernova rate. Finally the 

, morphology and kinematics of the circumnuclear molecular gas disks emerging from our models is 

■ similar to that revealed by recent near infrared observations using VLTI/Keck. 

\ Subject headings: galaxies: Seyfert - galaxies: starburst - ISM: structure - ISM: molecules - method: 

\Q . numerical 

o 
a^ 
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1. INTRODUCTION 

Optically thick molecular materials have been postulated to unify various observational properties of active galactic 
nuclei (AGNs), especially the two major categories of the AGN, namely type-1 and type-2 (Antonu cci 1993). It is 
^ often assumed that the obscuration is caused by dusty molecular gas ar ound a supermassive b lack hole (SMBH), 
■ whose geometry is schematically described as a doughnut-like 'torus' fe.g. lUrrv fc Padovanilll995l ). Its spatial extent 
and structures are deduced from comparison between observatior is in near /far infrared and theoretical spectra l energy 
distributions (SED) based on radiative transfer calculations (e.g. iFritz. Franceschini. fc Hatziminaoglou Il2006t lElitzud 
[2008). Detailed three-dimension al radiative transfer calculations suggest that the torus should not be uniform but 
clumpy (jSchartmann et al.ll2008[ ). In these model calculations, structures of the torus and its internal clumps are 
freely assumed, and then are tested in comparison with the observed SEDs. For example, the inner and outer edges 
of the torus are assumed typically as a sub-parsec and several tens of parcecs, respectively, and several hundred dense 
internal clumps are also assumed. It was then suggested that SEDs with the silicate absorption feature at 10 /im are 
mostly affected by the inner pc region of the torus (Schartmann et al. 2008). 

Recently the dusty gas in the circum nuclear region are directly observed in some nearby Seyfert galaxies using 
interferometric mid-infrared observations. The prot otypical type-2 S eyfert NGC 1068 harbors a warm (320 K) dust in a 
structure 3.4 pc in diameter and 2.1 pc in thickness (jJaffe et al.ll2004f ). VLTI/MIDI survey of nearby AG Ns with various 
types rece ntly revealed that in t he most objects AGN -heated dust is spatially compact (~ a few pc) (|Poncelet et aD 
[2OO6; Tris tram et al.l [2007t iMeisenheimcr et al.ll2007t 1 Tristram ct al., 2009,') . These observations showed a structure 
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of the hot dust and gas in a central few pc region, which may not necessarily encompass the whole structure of the 
molecular gas in the central region. Colder mo lecular gas p r obabl y extends further out, in i? ^ 10 pc to sub-kpc regions. 
Indeed using ^^CO (2 — 1) line observations, iHsieh et"an (|2008l ) revealed that the molecular gas in Seyfert 1 galaxy, 
NGC 1097, is concentrated in R ^llb pc, with an average density of the molecular hydrogen of riHa ~ 3 x 10"* cm~^. 
Central concentrations of high density gas around AGNs have been found via the observations o f molecular lines with 
high critical densities such as HCN J = 1 - (ricrit = 2.3 x 10^ cm'^) (e.g. iKohno et all 120031 ). 
If such quantities of circumnuclear molecular gas are present, then star formation is naturally ex pected. In fact, there 



are m any lines of evidence of ci rcum nucl ear sta rbursts in various types of AGNs arid quasars (|Ci d Fernande s et al 
20041: flmanishi & Wada 2004; D avies et al . 2007; Maiohno et aI1l2007tlRiffel et al.ll2007l:lWatabe eta l. 2008; Chen et al 



2009 ). iHicks et al.l (|2009l) recently revealed structures of the ISM in local AGNs traced by S{l)v = 1—0 line of molecular 
hydrogen at 2.1/xm using VLT/SINFONI. They found that the ISM whose scale is r ~ 30 pc forms a rotating disk 
with a high velocity dispersion, and correlated with the star formation rate. The average gas mass is estimated to 
~ IO^Mq, with a column density of TVh > 2 x 10^^ cm~^. They suggest that star formation occurs within a clumpy 
molecular gas phase. It has b e en observationally suggested that the starburst phenomena is linked with the nuclear 



activities ( Watabe et al . 2008; Chen ct al. '2009*), an d this connection was analytica lly studied base d on a clumpy , 



turbulent accretion disk (VoUmer & Bcckcrt 2003; VoUmer et al.|[200l iKawakatu fc Wada 2008; VoUmer et all 120081 : 
iCollin k Zaiml[200l . 

In order to construct physical models of the obscuring material around a SMBH, influenced also by concomitant 
star formation, three-dimensional radiative magneto-hydrodynamic calculations with a sufficiently high resolution and 
a fairly large dynamic range (e.g. from 100 Schwarzschild radii to 100 pc) are u ltimately ne cessary. How ever, until 
now, there are a few hydrodynamic studies i n the central tens pc (|Wada fc Nor man 2002; Wada fc Tomis aka 200^; 
lYamada et alll2007l: [Sctiartmann et al.ll2009D . IWada fc NormanI (|2OO0) (hereafter WN02) investigated evolution of the 
ISM in the inner 100 pc region around a 10* Mq SMBH using three-dimensional Euler-grid hydrodynamic simulations, 
taking into account self-gravity of the gas, radiative cooling and heating due to supernovae. They found that a clumpy 
torus-like structure is naturally reproduced on a scale of tens pc around a SMBH. The internal density, temperature 
are highly inhomogeneous, and the velocity field is turbulent, but the global geometry, i.e. a flared thick disk, is 
quasi-stable, implying that energy input from supernovae is balanced with turbulent dissipation. They suggested that 
AGNs could be obscured by the surrounding i nterstellar material. These results support observational findings that 
some AGNs are obscured by nuclear starbursts (jLevenson. Weaver, fc Heckmanll200lt iLevenson et al.l[2007t iBallantvni 
l2008l and references therein). However, the models in WN02 are not necessarily applicable to all circumnuclear regions 
with active star formation, because 1) the assumed supernova rate is extremely high (^ 10~* yr~^ pc^^), which would 
be appropriate only for intense s tarburst, like ultra-luminous infrared galaxies, but not to standard circum nuclear 
starbursts (e.g. iHicks et al.ir2009( ). and 2) the cold, dense gas in their simulations does not necessary represent the 
dusty molecular gas phase around an AGN, because H2 formation and its radiative destruction by far ultra-violet 
radiation (FUV) from massive stars and coUisional destruction are not included. As a result, the radiative cooling in 
WN02 is not consistent with the chemical abundances in the low temperature medium. 

Including the molecular gas in numerical models will eventually yield a much richer interface with current and future 
observations. This is because unlike the dust continuum, the effectively optically thin molecular line emission (even 
large line-center optical depths allow full view of turbulent media), is a direct probe of the velocity fields and thus of the 
gaseous disk dynamics, unobscured by the concomitant dust. The upcoming commissioning of ALMA will drastically 
enlarge such an interface by enabling the imaging of a whole suite of molecular lines probing different density and 
temperature regimes at high angular resolution. However from the results of WN02 it is still unclear how the various 
density regimes of the H2 gas would be distributed in the circumnuclear region, and what would influence its structures. 
We hereby extend the three-dimensional hydrodynamic simulations of WN02 to explicitly track the evolution of the 
H2 phase in the central 64^ x 32 pc region, and its interplay with the HI phase, with a 0.125 pc resolution (cf. 0.25 
pc in WN02). We also vary the strength of the uniform FUV field and the supernova rate (SNR) in order to examine 
their effect on the structures of molecular gas. 

2. NUMERICAL METHOD AND MODELS 
2.1. Numerical Methods 

We use a numerical scheme based o n Eule r ian h ydrodynamics with a uniform grid, which is the same as those 
described in IWada fc NormanI (j2001[ ): IWadal (j2001l ). Here we briefly summarize them. We solve the following 
equations numerically to simulate the three-dimensional evolution of a rotating gas disk in a flxed spherical gravitational 
potential. 

dp/dt + V ■{pv)^Q, (1) 

dv/dt + {V ■ V)V + Vp//?= -V$ext - V$BH - V$sg, (2) 

dE/dt + V • [{pE+p)v]/p^Tvy{Go) + FsN - pA(Tg, /h,, Go), (3) 
V2$,g = 4^Gp, (4) 

where, the speciflc total energy E = {v]"^ /2 + p/ {j — l)p, with 7 = 5/3. We assume a time-independent external 
potential $cxt = -(27/4)i/2[„2/(-^2 _^ ^2^1/2 _^ yljir^ + alYl\ where ai = 100 pc, a-i = 2.5 kpc, vi = 147 km s-\ 
V2 = 147 km s-i, and $bh = -GMBu/{r'^ + h'^f'^, where Mbh = 1-3 x lO^Afg and & = 1 pc. See Fig. [Jfor the 
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rotation curve based on the external potentials $cxt and <i>BH- Observationally it is hard to determine exact rotation 
curves in the central 100 pc region of external galaxies, the adopted mass distribution with a core radius of 100 pc i s 
roughly consistent with the rotation curves around nearby Seyfert nuclei observed by VLTI/Keck (|Hicks et al.ll2009l ). 
CO (1-0) observations of the central region of nearby spi ral galaxies also s how that their rotation curves are steep in 
the central part, suggesting central massive c omponent s (Sofue et al.lll99 9). This could correspond to the stellar cores 
found by HST/NICMOS in nearby galaxies (jCarollo et al...200^ . The potential caused by the BH is smoothed inside 
r '--^ 1 pc by introducing a 'softening' parameter 6, in order to avoid too small time steps around the BH. This softening 
does not change the rotation curve in the range shown in Fig. [U but we should not rely on the gas dynamics and 
structures at r < 1 pc. In the central grid cells at r < 0.25 pc, physical quantities stay constant to represent a sink 
cell. 

The hydrodynamic part of the basic equation s is solved by AUSM (Advection Upstream Splitting 
Method) ((Liou fc SteffenlTl99a I Wada fc Nor"^^[200l . We use 512^ x 256 grid points in high resolution models 
and 256^ x 128 grid points in low resolution models to investigate long term evolution with different parameters^. 
Cartesian grid points cover a 64^ x 32 pc'^ region around the galactic center (i.e. the spatial resolution is 0.125 pc in 
the high-resolution runs and 0.25 pc for the low resolution runs). The Poission equation, eqn. (|4|) is solved to calculate 
the self-gravity of the gas using the Fast Fourier Transform and the convolution method, where 1024^ x 512 grid points 
and a periodic Green's function is used to calculate the potential of an isolated system (jHocknev fc Eastwoodlll98 11. 

An essential progress in the present calculations compared to WN02 is that we now solve non-equilibrium chemistry 
of hydrogen molecules along with the hydrodynamics. Formation of H2 on dust and its radiative and coUisional 
destruction are explicitly tracked, therefore we can deduce the distribution of H2 in the central tens pc region. See 
details in Appendix A. 

In the energy equation (eg. (3)) , we use a cooling function based on a radiative transfer model of photo-dissociation 
regions (|Meiierink fc SDaansl[2005h . K{Tg, fu^,Go) (20 K<Tg< 10* K) assuming solar metallicity, which is a function 
of the molecular gas fraction /h2 and intensity of FUV, Go (see Appendix B for the detail). The radiative cooling rate 
below ~ 10^ K is self-consistently modified depending on /h2, which is determined in each grid cell. We adopt heating 
due to the photoelectric effect, Fuv and energy feedback from type-II SNe (see below), Fsn- We assume a uniform 
UV radiation field, whose strength is an important parameter in the model, represented by the Habing unit (Go), the 
incident FUV field normalized to the local interstellar value. 

2.2. Initial conditions and Model parameters 

The initial condition is an axisymmetric and rotationally supported thin disk with a uniform density profile (thickness 
is 2.5 pc) and a total gas mass of Mg = 6 x lO^M©. Since we allow outflows from the computational boundaries, the 
total gas mass decreases during the evolution, and it settles to ~ 5 x lO^M© at a quasi-cquilibrium state. Random 
density fluctuations, which are less than 1 % of the unperturbed values, are added to the initial gas disk. A uniform 
disk with p = 1.22 x IQ^Mq pc~^ at t = is evolved for t = 1.4 Myr (the rotational period at i? = 10 pc is about 0.2 
Myr) using 256^ x 128 grid points (i.e. spatial resolution is 0.25 pc), then it continues to higher resolution runs with 
512^ X 256 grid points. The system settles into a quasi-equilibrium state, where the total molecular gas mass is nearly 
constant, after t ^ 3.5 Myr (see §3.2 and Fig. (5]). 

Since the current integration time is short (~ 5 Myr) due to a limitation of our computational resources, instead 
of tracking the whole life of massive stars, here we assume that supernova (SN) explosions occur at random positions 
with a constant rate in the region confined by r < 26 pc and jz] < 2 pc. Here we effectively assumed that massive, 
cold molecular gases from which stars formed are not distributed spherically, but they are rather concentrated in a 
disk with a small scale height. Therefore we expect that supernova explosions mostly occur near the disk plane. 

The average SN rate (SNR) is one of the free parameters in the present study. Observationally there is a large 
ambiguity in the star forming activities around AGNs. the star formation rate (SFR) around AGNs was suggested to 
several tens M0yr~^ kpc~^ (jPavies et al.ll2007l) . For the Salpeter initial mass function (IMF) and stellar mass rage 
of 0.1-125 Mq, type-II supernova rate is ~ 0.007 SFR, and a corresponding supernova rate per unit area would be 
^ 10~^ yr~^ pc~^. For the circumnuclear disk of r = 26 pc in the present model, this rate suggest that ~ 3 x 10~"^yr~^. 
On th e other hand, using the scaling relation on the star formation rate in galactic disks and in starbursts ([Kennicuttl 
Il998| ). the initial gas density in the present model yields 4 x 10~* yr~^- We here change the supernova rate from 5.4 
to 540 xlO~^ yr. For each type-II supernova explosion, the energy of 10^^ ergs as thermal energy and gas mass 8 
Mq are instantaneously injected into a single cell. Note that since we set a maximum temperature of 10^ K, not all 
the injected thermal energy is used to create blast waves. Three dimensional evolution of blast waves caused by SNe 
in an inhomogeneous and non-stationary medium with global rotation are followed explicitly, taking into account the 
radiative cooling, which makes the evolution of the supernova remnants depending on the local gas density distribution 
around the SNe. 

We here assume a spatially uniform FUV field with Go = 10 and 100. For a p hoto-dissociation region (FDR) 
surrounding an HII region, Go ~ 2 X 102(L/10'*LQ)(r/l pc)-^ (e.g. iTielensI [20051 ) . For the star formation rate 
lO^^Af0yr~^ pc~^, which is a typical rate in starburst galaxies, the luminosity of massive stars formed for 10^ yr that 
contribute to the FUV field is roughly lO'' L0pc~^. If we assume that the massive stars are embedded in the uniform 
gas disk with n ^ 10^ cm~'^, contribution from local (< a few pc) stars is important. Then the average local FUV field 

^ In the high resolution models, a calculation for 4.5 Myr requires about 180 CPU hours using 16 CPUs of NEC SX-9 whose peak 
performance is about 1.6 TFlops. 
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Go 


A fncl 


SNR [lO"-'' yr"^] 


HlOa 


10 


0.125 


5.4 


LlOa 


10 


0.25 


5.4 


HlOOa 


100 


0.125 


5.4 


LlOOa 


100 


0.25 


5.4 


LlOOb 


100 


0.25 


54.0 


LlOOb* 


100 


0.25 


54.0 


LlOOc 


100 


0.25 


540.0 



TABLE 1 

Model parameters. "H/L" represents "High/Low" resolution (A, size of a numerical grid cell), "100/10" represents 

INTENSITY OF FAR ULTRA-VIOLET RADIATION IN THE HaBING UNIT (Go), AND SUFFIX 'A-C' REPRESENT THE AVERAGE SUPERNOVA RATE 

(SNR). LIOOb* is the same as LIOOb, but energy from supernovae is injected in a larger scale height, i.e. \z\ < 10 pc. 

in a grid cell of 0.1 pc would be Go ^ 200. On the other hand, an AGN would have a negligible effect on the FUV 
field in the most regions, because the average column density through the initial thin disk is > 10^^ cm~^ for r > 2 pc. 
However, as we will see in §3, the star forming disk in a quasi-equilibrium state is highly inhomogeneous, especially in 
the outer region > 5 pc). Therefore, Go could be several 10'^ for the low column density regions {N < 10^^ cm~^) 
at r ~ 10 pc from a typical AGN with luminosity L = lO''^^ erg s~^. In reality, both the density field and distribution 
of massive stars are not uniform, and Go should have a large dispersion with some radial dependence. Effect of the 
non-uniform radiation field will be an important subject to solve by three-dimensional radiative transfer calculations 
in the future. 

Beside the FUV, the X- ray radiation from the AGN would be impo rtant for the thermal and chemical structures of 
the ISM around an AGN (jMalonev et al.lll996t FMeiierink et al.ll2007| ). This wiU be discussed in §5. 

Model parameters are summarized in Table 1. For simplicity, we here assume that the FUV intensity and SNR are 
independent free parameters. In reality those quantities are tied to the ambient star formation rate density, which 
in turn is regulated by the H2 gas density since stars always form out of the H2 gas. It is therefore expected that 
FUV, SNR and H2 gas mass fraction evolve in a highly coupled manner, leading perhaps towards a self-regulating, 
quasi-stationary state. This important issue can only be explored if a SF-regulating H2 gas distribution is included in 
the ISM-|-stars model, and we intend to do so in a future paper. 
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Fig. 1. — Assumed rotation curve (thick solid line) and contribution from each component of the external potential. Thick dashed line 
(ai = 100 pc, Di = 147 km s~^) and thin line (02 = 2500 pc,t'2 = 147 km s~^) are from 3>ext, and thick dotted line is from 'I'bh with 
A/bh = 1.3 X IO^'Mq. 

3. RESULTS 

3.1. Structures of a fiducial model 

In a quasi-steady state {t > 3.5 Myr), as reported by WN02, the gas forms a highly inhomogeneous, flared disk, 
as seen in Fig. [2l which shows gas density, temperature and molecular hydrogen density in model HlOOa at i = 4.38 
Myr. The temperature map shows that cold {Tg < 100 K) gas is mainly distributed in the high density regions. In 
the central funnel-like cavity, the temperature of the gas is hot {Tg > 10^ K). Hot gases are also patchily distributed 
in the cold, flared disk. Typical size of these hot cavities is a few pc. A large fraction of the volume is occupied by 
warm gas {Tg ~ 8000 K). As expected distribution of H2 roughly follows the cold, dense gas, and therefore it forms a 
high density circum nuclear disk whose radius is about 5 pc, surrounded by a porous torus which extends to ^ 5 — 10 
pc above the disk plane. 

Figure [3] shows velocity fields of two components {vy and Vz) at the same snapshot of Fig. [S] The torus shows a global 
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Fig. 2. — Cross sections (a) gas density [Af0pc '^], (b) temperature [K] and (c) H2 density [AfQpc ^] on x-y and x-z planes of model 
HlOa {t = 4.38 Myr). 

rotation, but there are also large internal random mot ions. This comp licated velocity field is naturally expected in an 
inhomogcncous, self-gravitating disk (Wada, Mcurer, fc Normaiill2002l ). but the vertical motions are mainly enhanced 
by energy input from supernovae (see also Fig. (H). The vertical velocity field also shows a bipolar outflow in the 
central funnel, which is mostly warm and hot gases as seen in Fig. ^p. 

In Figs. and|4)D, we plot total column density {Ng) and H2 column density (N-a^) toward the galactic center as 
a function of the viewing angle in model HlOa. It clearly shows that the total gas and H2 column densities are both 
largest at 0^ ~ (i.e. edge-on), and the average column density decreases on average towards the pole-on view as 
expected. We should note, however that the scatter around the average is two orders of magnitude or more for any 
viewing angles, which is comparable to the change of the average value between 9y = 0° and ±90°. Figure [Hd shows 
II2 is distributed similary, but it is more concentrated near the disk plane. The scatter of Nn^ is more significant than 
that in A^g especially for 9 > 50 deg, reflecting that II2 gas is highly inhomogeneous and more sparsely distributed in 
larger latitudes. 

3.2. Dependence of Model Parameters 

Figure 5 shows density distributions of models with two resolutions, HlOOa (A = 0.125 pc) and LlOOa (A — 0.25 
pc) at t = 4.02 Myr. The inhomogeneous density fields and global morphology of the thick disk in the two models are 
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Fig. 3. — Same as Fig. [2] but for velocity fields of the gas: (a) Velocity for ?/— direction, and (b) velocity for 2:— direction. 

qualitatively similar, but clumpy structures are finer in the high resolution model, HlOOa than in LlOOa. In Fig. [6] 
evolution of molecular gas fractions, = -^Ha/Aftotai in five models are shown. The molecular gas fraction in HlOa 
rapidly decreases from — 0.5 to ~ 0.32 in the first Myr, then it increases and slightly oscillates around /hj — 0.4. 
The model with ten times stronger FUV (HlOOa) shows a similar time evolution, with f^^ only about 10% smaller 
than in model HlOa. In order to see a longer term evolution, we restart calculations from a snapshot of the high 
resolution data at i = 2.5 Myr with a lower resolution (i.e. 0.25 pc)^. The low resolution run (LlOa) shows rapid 
decrease of /h2 in the first ~ 0.2 Myr, reverses to an increase at i ~ 2.7 Myr, and then settles around — 0.4 after 
t ^ 3.5 Myr in LlOa as in HlOa. Molecular fraction of LlOOa in the quasi-equilibrium state is /hj — 0.32. These 
results suggest that for a given supernova rate and the FUV intensity there is an equilibrium H2 fraction which does 
not significantly depend neither on the numerical resolution nor on the initial state. For stronger FUV fields, H2 is less 
abundant, but this is not very sensitive to Gq because of the effective H2 self-shielding against far-UV dissociation. 
In the low resolution runs, we change energy input rate due to the supernova rate by a factor of 100 (see Table 1). 
As seen in Fig. [6l LlOOa and LlOOc show similar evolution, but the molecular fraction is only slightly smaller in the 
model with 100 times larger supernova rate. 

Figure [7] shows column density distribution of II2 in models HlOa and HlOOa at t = 3.47 Myr. The dependence 
of iYjjj on the viewing angles are similar in the two models, but in HlOOa, H2 more rapidly decreases toward high 
latitudes. This is naturally expected because H2 is mainly formed in the dense regions, which are concentrated near 
the disk plane. Strong FUV can easily dissociate H2 in the relatively diffuse regions at high latitudes. However, those 
differences in the two models are not significant. Even in the strong FUV model, H2 is inhomogeneously distributed 
in the central tens pc. 

We compare vertical velocity dispersion of the gas below Tg = 2000 K in three models (LlOOa, LlOOb, and LlOOc) 
with different supernova rates in Figure [5] It shows that velocity dispersion tends to increase with higher supernova 
rate. However its dependence in the torus region is less than factor of two when the supernova rate is 100 times larger. 
The effect is more significant in the inner region, where there is a bipolar outfiow from the central funnel (Fig. [3]). 

4. DISCUSSION 

4.1. Comparison with Observations 

iHicks et al.l ()2009l ) recently revealed structures of the ISM in local AGNs traced bya5(l)i' = l — line of molecular 
hydrogen at 2.1/im using the near-infrared field spectrograph SINFONI at ESO Very Large Telescope and OSIRIS 

^ The computational time in the high resolution run is about 16 times longer than that in the low resolution runs. 
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Fig. 4. — (a) Column density of the gas, A^g and column density of H2 , Nn^ as a function of the viewing angle in model HlOa at t = 4.38 
Myr. The solid line represent an azimuthally averaged column density. 

at Keck Observatory. They found that ISM with a radius of r ~ 30 pc forms a rotating disk with a high vefocity 
dispersion (several tens km s^^), and suggested that star formation is concomitant with clumpy molecular gas. The 
average gas mass is estimated to ^ IO^Mq, assuming a 10% mass fraction to the dynamical mass. The size, mass and 
dynamics of these nuclear molecular gas disk around AGNs are very similar to what we found here. The turbulent-like 
velocity field revealed by their observations qualitatively resembles to those seen in our models (Fig. [3]). They also 
suggest a positive correlation between the velocity dispersion of the hot molecular gas and the star formation rate, 
which is consistent with our result (Fig. [5]) , and it is reasonable if turbulent velocity is energized by supernovae and 
balanced by gas dissipation (WN02). Although 1/ = 1 — 05(1) is emitted from warm gas, it is expected from our 
si mulations that the cold molecular gas have also similar random velocity field. 

'Hs ieh et al.l ()2008f ) observed a Seyfert 1 galaxy, NGC 1097 using ^'^CO{ J = 2 - 1) by the SubmiUimeter Array at 
angular resolution of 4.1" x 3.1" (290 pc x 220 pc), and deduced a total mass of molecular hydrogen of Mn^ — 
6.5 X 10''' Mq, with column densities A^Ha.co — 5.9 x 10^^ cm^^, and A''h,co — 1-2 x 10^^ cm^^. They suggest a 
clumpy disk configuration to explain inconsistency between the large column density suggested by their observations 
and two orders of mag nitude smaller column density (A'h.x — 1-3 x 10^^ cm~^) suggested by X-ray observations 
([Terashima et al.l |2002|) . If the nucleus of this galaxy surrounded by a torus consisted of high density, small gas 
clumps, it would be possible to observe the broad line region through gaps between the clumps. On the other hand, 
the H2 column density represents an average one in the large beam size (~ 250 pc), which should be much larger than 
than the value suggested by X-ray observations. Could this be the case in the clumpy H2 in our results? The ratio 
■^H2,Co/-^H,x is — 45 in NGC 1097. Among our models, model HlOa, in which the moderate star forming activity 
is assumed (Go = 10 and SNR= 5.4 x 10~^ yr~^), would be appropriate for comparison. If one supposes that the 
H2 column density inferred from the CO line observations represents an average column density of all the molecular 
gas present in the central several tens parsecs {{NH2)), and column density of H suggested by X-ray reflects the amount 
of material along the line-of-sight (its viewing angle is 6^) for the nucleus {Nb_{Ov))- In Fig. [9l we plot (A'h2)/-^h in 
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Fig. 5. — Density distributions on x-y and x-z pianos (a) in a high-resolution model HlOOa and (b) in a low resolution model LlOOa at 
t = 4.02 Myr. 
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Fig. 6. — Evolution of molecular gas fraction /h2 in high resolution models (HlOa and HlOOa) and low resolution models (LlOa, LlOOa 
and LlOOc) recalculating from the high-resolution runs at t = 2.5 Myr. 

model HlOa at i = 4.38 Myr. It is clear that we have smaller chances to observe small iVn if the line of sight is closer 
to the edge-on. (A''h2)/^h ^ 10 can be achieved for 50 deg or larger. One should note again that the scatter 

of the ratio is significant for any given viewing angles, and it also depends on the structure of molecular gas around 
the nucleus. It is therefore hard to conclude at this moment on structures of the circumnuclear region of NGC 1097. 
Future observations of AGNs using CO or other lines with much higher resolutions, at least with a few pc beam, by 
for example ALMA are necessary. 
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Fig. 7. — Column density of H2 as a function of a viewing angle in HlOa (Left: Go = 10) and HlOOa (Right: Go = 100) at t = 3.47 Myr. 
The solid line represent an average column density. 
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Fig. 8. — Azimuthally averaged velocity dispersion for z-direction as a function of radius in three models at t = 5.32 Myr. Supernova 
rates are 5.4x10"^ yr"! (LlOOa), 5.4x10"" yr"! (LlOOb), and 5.4x10"^ yr"! (LlOOc). 
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Fig. 9. — (AfH2)/^H as a function of the viewing angle for model rd at i = 3.79 Myr. The horizontal line shows the ratio suggested by 
C O and X-ray observations in NGC 1097. 

iJafFe et all p004) observed the central region of NGC 1068 by the Very Large Telescope Interferometer (VLTI), 
and suggested a dusty "torus" of thickness 2.1 pc, diameter 3.4 pc, and temperature of Tg = 320 K. This could 
correspond to the central, relatively smooth molecular dense disk typically seen in our models (Fig. [2]). Due to the 
strong gravity and shear caused by the central BH, it is natural that molecular gas forms a thin, smooth disk in the 
central few pc. However, if the observe d central d isk of NGC 1068 is geometrical ly thick, a stron g heating source, e.g. 
infrared radiation would be necessary ( [Pier fc Kr olik 1993; Krohk 2007). Tristr am etldl ()2007D suggested a clumpy 
torus plus a central thick disk model from mid-IR observations using VLTI of Circinus, which has a Seyfert 2 nucleus. 
Although the size is smaller (~ 2 pc), the geometry and structures suggested by their high resolution observations (see 
a schematic picture in their Fig. 9) is quite similar to the models presented here. 

4.2. Effects of X-rays 

The strong X-ray continuum emanating from th e AGN could affect the chemical state of the ISM in the central 
region (iMalo ney et al.l 119961 : iMeijerink et al.l l2007l ) . It is known that the chemical state of the X-ray Dissociated 
Region (XDR) is mainly determined by Hy^/n^ where iJx is a X-ray heating deposition rate and n is the number 
density of the gas (Malone v et al.|[l996l ). Although we do not include the effect of X-ray from the nucleus, here we 
estimate its potential effect using one of our models. In Fig. 1101 we plot a fraction of II2 and temperature as a 
function of Hyi/n (erg cm'^ s~^) at randomly selected points in the computational box. Here we assume that the X-ray 
luminosity of the ACN is Lx = 10'*'* ergs s~*. Although there is a large scatter, for \og{Hx/n) ~ —24 or smaller the 
molecule fraction increases significantly and the temperature in most regions is less than 100 K. At high gas densities 
(~ 10^ cm~^), large (> 0.1) H 2 fractions can be sustained for the range of Hx/n values, limiting the impact of XDR 
effects (|Meiierin k et al.l 120071 ). This is partly because the free electrons associated with ionization processes help to 
form H2 through the H~ route and also because H2 formation scales with density squared. Still, for Hx/n significantly 
larger than 10~^^, gas is mostly atomic as expected (Maloney et al. 1996). 

These results suggest that XDR chemistry may change distribution of H2 around an AGN, if we explicitly include 
X-ray from the nucleus in the following sense. The H2 abundance is indeed robust in a clumpy mediu m, but the 
XDR gas temperature should rise relative to the PDR case by a factor of ^ 5 for log Hx/n > —26 (Mei jerink et al.l 
l2007f ). After all. X-ray ionization heating with for example Hx/n ~ 10~^^ erg s~*, is more efficient than photo-electric 
emission by dust grains, therefore H2 gas would have higher temperatures. This leads to stro nger emission in the 
pure rotational H2 hues, like e.g. S{0), S{1) at 28/xm, 17^m. Similar considerations hold for CO (jSpaans fc MeiierinS 
[2OOI . Shocks can heat as well, therefore observational data on line profile kinematics would be useful in the future. 

4.3. A self-regulated H2 and the star formation rate 

In the present results it is notable that although the assumed FUV and supernova rate differ by orders of magnitude 
among the models, total molecular gas fraction settle between /hj — 0.3 — 0.4. There is a weak 'negative' feedback 
of FUV and supernovae on the star formation in the central region. This suggests that the molecular gas disk could 
survive around AGNs associated with a starburst. Since the negative feedback of FUV and supernova rate on H2 mass 
is 'weak', it is also implied that if there is an episodic large mass inflow into the circum nuclear region, we expect a burst 
of star formation, because star formation rate de pends on H2 den sity. This could also affect the mass accretion rate 
toward the AGN via turbulent viscosity (Kawakat u fc Wadall2008f l. To explore these mechanisms, the star formation, 
the radiative feedback from t he AG N /stars affecting H2 (the star formation fuel), and a coupled SNR should be 
self-consistently followed (e.g. IWadali2008 ) for at least several tens Myr. These effects, along with incorporating the 
H2 -tracing CO molecule, will be investigated in a future paper. 
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Fig. 10. — (a) H2 fraction as a function of ratios of the X-ray heating deposition rate and gas density, Hx/nlerg cm^ for randomly 
selected 500 points in the computational box in model HlOa at t = 4.3 Myr. (b) Same as (a), but for temperature. Note that H2 fraction 
is zero in many selected points, because the hot and warm gas occupies a large volume, therefore the number of data points in the plot (a) 
is smaller than in (b). 

4.4. Positions of supernovae 

Finally, effect of position of supernovae on the gas density distribution is compared between LlOOb and LlOOb* (Fig. 
Ilip . As explained in §2.2, we assumed that SN explosions occur at random positions with a constant rate in the region 
confined by \z\ < 2 pc. This is justified by the results that H2 is more concentrated in the disk plane. However, orbits 
of massive stars formed near the disk plane could have large inclination angles, as a result, energy from supernovae 
could be released at higher latitudes. In model LlOOb*, supernovae are assumed to be exploded in \z\ < 10 pc (Fig. 
[TTb). but we do not see significant differences in the distribution of gas around AGNs. 




Fig. 11. — Density distributions in models with different positions of supernovae, \z\ < 2 pc (LlOOb) and \z\ < 10 pc (LlOOb*). 



5. CONCLUSIONS 

We present new high resolution numerical simulations of the interstellar medium (ISM) in a central i? < 32 pc region 
around a supermassive black hole (1.3 x IO^Mq) at a galactic center. Three-dimensional hydrodynamic modeling 
of the ISM (Wada & Norman 2002) along with a concomitant starburst now includes tracking of the formation of 
molecular hydrogen (H2 ) out of the neutral hydrogen phase as a function of the evolving ambient ISM conditions 
with a finer spatial resolution (0.125 pc). Radiative cooling rate is sclf-consistently changed depending on abundance 
of II2 in each grid cell. In a quasi equilibrium state, mass fraction of H2 is about 0.4 (total II2 mass is ~ 1.5 x 10^ Af©) 
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of the total gas mass for the uniform far UV (FUV) with Gq — 10 in Habing unit. As shown in the previous model, 
the gas forms an inhomogeneous disk, whose scale-height becomes larger in the outer region. H2 forms a thin disk in 
the inner ~ 5 pc, and it is clumpy and swelled up toward /i ~ 10 pc. The velocity field of the disk is highly turbulent 
in the torus region. Outflow of the gas forms a hot (Tg ~ 10^ K) funnel in the inner region. These structures and 
dynamics are energized by supernovae with average supernova rate (SNR) of ~ 5.4 x 10~^yr~^. Gas column densities 
toward the nucleus larger than 10^^ cm~^ are found if the viewing angle is smaller than 0„ ~ 50° from the edge-on. 
However, the column densities are distributed over almost two orders of magnitude around the average for any given 
viewing angle due to the clumpy nature of the torus. The column density of H2 shows a similar distribution, and 
A'h2 10^^ cm~^ is achieved for < 30°. For stronger FUV (Gq = 100), the total H2 mass in an equilibrium is only 
a little smaller (~ 0.35) and the molecular gas is slightly more concentrated in a plane, while other properties of the 
ISM are also not very sensitive to the FUV intensity and supernova rate. Vertical velocity dispersion of the gas in the 
torus (~ 20 km s^^ on average) is enhanced only by a factor of ^ 1.5 for increasing the supernova rate by two orders 
of magnitude. 

The structure and dynamics of these nuclear molecular disk presented here are very similar to molecular gas disk in 
various types of n earby AGNs recently revealed by ro-vibration line of H2 using VLTI/SINFONI and Keck/OSIRIS 
(jHicks et al.ll2009t ). The inconsistency on column densi ty toward t he center in NGC 1097 between the large column 
density suggested by CO observations (jHsieh et al.ll2008l ) and X-ray (jTerashima et al.l[2002[ ) could be attributed to the 
inhomogeneous structures emerging from the hereby presented disk models. 

Finally, in the present results it is notable that although the assumed FUV and supernova rate may differ by orders 
of magnitude among the models, the total molecular gas fraction settles between — 0.3 — 0.4. This weak 'negative' 
feedback from FUV and supernovae on the star formation in the central region, suggests that substantial molecular 
gas disks can survive around AGNs with strong ongoing star formation. 
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APPENDIX 

NUMERICAL RECIPE OF H2 FORMATION AND DESTRUCTION 

Abundance of hydrogen molecule of each hydrodynamic grid-cell (0.125 or 0.25 pc) at a given time step is determined 
by the following time-dependent approach coupled with hydrodynamic equations. In this approach H2 and HI are 
assumed fully mixed within "cells" of typical radius Rc whose choice is motivated by the ISM physics, while it remains 
unresolved by the simulation during the course of the thermodynamic and dynamic evolution of the gas. This allows a 
smooth treatment of the III-II2 mass exchange in both the Cold Neutral Medium (CNM) and Warm Neutral Medium 
(WNM) phases, with all factors now incorporated in a single equation, namely 



dn2 (r) 
dt 



Ri{T^)n{rf 
2i?f(Tk)n(r) 



47r 



fs{N2)e-^^'''^dn 



4ir 



-7i(Tk) [n{r) - 2n2(f)] n^ir^ - 72(Tk)n2(f)2. (Al) 
^ is the unshielded FUV-induced H2 dissociation rate, Go is the FUV radiation field in 



-11 , 

^ ergs cm~^ t is the FUV optical depth, n{r) = ni(r) -I- 2n2{r) [rii: neutral hydrogen 



The factor /cd = 4 x 10 
Habing units (1.6 x 10" 

density, 712: II2 density), and 71 (Tk), 72 (7k) are t he H-H2 and H2 -H 2 colhsional destruction rates of II2 for gas 
temperature Tk that can be fitted using data from iMartin et all (|1998l ). The function i?f(rk) is the II2 formation 
factor. We set i?f = 3.5^* x 10-i^(rk/100 KY/^Sb_{T\,) cm^ s'^ (IJuralll975t fPelupessv et al.l[2006h for Tk < 500 K, 
otherwise Rf ~ 0. The factor /i represents obsevational uncertainties in the H2 formation rate, with the canonical value 
(fjura 1975) corresponding to /i = 1, but values up to /i = 5 are possible (sec Papadopoulos et al. 20 0^). Here we adopt 
ji = 5, since this seems more probable form ISO studies of Galactic Photodissociation regions (jPapadopoulos et al.l 
I2OO2I and references therein). The function SuiTi^) quantifies the sticking probability of an H atom on a grain, which we 
set to unity for most of the modeling done in this work (the H2 formation probability after sticking is considered unity 
throughout). We only test the case of = (1 Tk/To)" where To = 100 if which was obtained bv lBuch fc Zhangi 
(| 199 If ), and is thought to be good for dust grains in general, as long as they are covered by a few layers of weakly 
bonded material. 

The /s(iV2) function denotes the H2 self-shielding factor, with A^2 being the H2 column density, and 



/s = l, for N2<Nn^ lO^^ cm-^, 




where fc=3/4 fsee lDraine fc Bertoldil (|1996D ). 

If we integrate eq. (Al) over the volume AVc of the cell, while assuming negligible dust absorption within it 
(r(ri) ^ 0) we obtain, 



dn2 

~dr 



2Rf{Ti,)n + Gokd 



No 

712 i?c 



K{Rc 



n2 



7l(7k) ("- - 2n2) 77.2 - 72(Tk)7ll, 

where the factor K(Rc) {Rc is a cloud radius) is 



(A3) 




-fc+2 



] dx, where x = r/Rc, 



K{Rc) = TrTT7 I I {ttt^] dnd\- (A4) 

(A5) 

(3-fc)(4-2fc) [ 4-fc ' ^^^^ 

and for k — 3/4 results to K{Rc) = 0.96. Unlike [Pelupessv et al.l (120061 ). the description of the HI-H2 gas mass 
exchange here does not involve the density-size power law. It utilizes only the much better established line-width-size 
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scaling law found for molecular clouds (jHeverfc Brunij [2004D . and used here to define the "cell" size (A. 7) within 
which HI and H2 are assumed fully mixed. Aside from the much smoother implementation of the H2 formation and 
destruction mechanisms that this allows, it yields a treatment that remains suitable for a gas phase that may not 
follow both scaling relations (i.e. diffuse and non self-gravitating gas). 

The assumption of a typical HI-H2 mixing sub-grid "cell" where H2 formation and dissociation takes places in full 
concomitance with HI greatly simplifies the problem of tracking the HI<^H2 gas phase exchange by reducing a non- 
local radiative transfer problem to a local one. Thus the H2 fractions estimated are likely to be lower limits since 
H2 self-shielding from FUV radiation occurring between distant "cells" (or dust absorption of the FUV radiation over 
scales larger than 2Rc) can only result to higher H2 fractions. 

The local approximation made here in order to simplify a formidable radiative transfer problem in turbulent media 
still leaves the choice of the HI-H2 mixing scale i?c open. A natural choice of Rc for our three dimensional hydrody- 
namic-|-thermodynamic ISM models would be Rc ~ 1/2 Ip, where Ip is the scale where the turbulent ISM pressure 
begins to dominate the thermal pressure, or equivalently the scale where the cloud size i?c-liiiewidth {AV) relation 
breaks down and AV{lp) AVth- The corresponding radius is then given by: 



Rc = 0.5 



AVt 



th 



1/9 



PC with AVth =i 0.8 (Tk/lOO K)i/2 km s-\ 



(A7) 



where AV{1) — AVo{l/pcY i s the size-linewidth relation found for the cool ISM. Several observational and theoretical 
studies favor q ^ 1/2 (e.g. 'LarsonTosT; 'Elmegrcen '1989'; He verfc BruntI [2004f) . The normalization factor AVq = 
1.2 km s^^ for both CNM and WNM (Wolfirc ct al. 2003). Below such scales the thermal velocity fields take over, 
setting a natural resolution limit for any hydrodynamic ISM models tracking macroscopic gas motions. Moreover 
tracking the thermal state of the gas as a functio n of time and position allows Rc t o adjust accordingly (cf. eq. IA7|) . 
as expected in the true multiphase medium (e.g. IChieze &: Pineau Pes Foretslll987D. Fi nally, by avoiding linking Rc 
to the numerical resolution of our hvdrodvnamic models (e.g. iGlover &: Mac Lowll2007f l keeps the computed H2 gas 
mass fr actio n independent of it. We set a maximum value of Rc corresponding to the scale length of Tk ~ 500 K gas 
(cf. eq. IA7[) . This temperature limit broadly marks the CNM^ WNM thermodynamic transition zone beyond which 
purely thermal motions fully overtake the macroscopic ones over the scales resolved here, while little H2 content is 
expected (i?/ ~ for Tk > 500 K while colhsional destruction of H2 becomes very efficient). 

In the case of CNM-to-WNM or the WNM gas phase Rc can be substantial, and the dust extinction term can be 
important. Then we use 



K{Rc) = 



Rc 



which after some simplifications results to, 



1+X 



^l-feg-(na-Fuv-Rc)t^^ 



dx, where x ~ r/ Rc 



(A8) 



(A9) 



The latter expression encompasses both dust extinction and H2 self-shielding and is the correct one for use in eq. (A3). 
For apuv — 0, K{Rc) becomes identical to that given in eq. (A6). In the case of pure dust shielding (i.e. k — 0), and 
after some calculations it becomes 



K{Rc) 



{2a - 3) + (2a2 + 4a -f 3)6" 



, a — n apvyRc 



(AlO) 



COOLING FUNCTION BASED ON PDR MODELS 

The cooling curves in temperature (Tg), H2 abundance and radiation field strength Go are derived from a grid of 
photo-dissociation region (PDR) models that have a range in irradiation (Go = 0—10^ values) and density 1 — 10® cm^^. 
The radiation field extends beyond 13.6 eV acco rding to a starburst9 9 spect rum for a Salpeter IMF, i.e., the HII region 
is computed as part o f the PDR. The code of iMeiierink fc SpaansI (|2005f ). with the latest chemical rates and grain 
surface reactions as in ICazaux &: Spaani ()2009f ) is adopted. 

The PDR models, for a given Go and density, are evaluated at the particular H2 abundances and temperatures 
that they enjoy with depth to get a parameterization of the cooling rate in Go, H2 abundance and Tg. In general, 
the chemical and thermal balance in PDR models depends on the ambient density. As it turned out, a cooling curve 
extraction (density squared scaling) was possible because many cooling lines have critical densities well in excess of 
10^ ciR-^. 

However, optical depth effects can suppress cooling and drive level populations to LTE. Therefore, information on 
the typical column over which coolants exist chemically in a PDR are used to correct for line trapping, as follows. A 
local turbulent velocity dispersion (Gaussian a) of dV = 3 km s~^ is adopted and some turbulent coherence length 
L (between 3 and 0.3 pc). This yields a formal velocity gradient dV/L for scales larger than L. If the opacity is 
concentrated on small physical scales, then radiation is trapped. Still, different parts of the PDR cloud do not see each 
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other on large physical scales. In effect, one has an effective Sobolev approximation or not depending on the chemical 
stratification. We adopt L = 0.3 pc in these sim ulations to mimic th e effec ts of a turbulently active nuclear region. 
In this, the multi-zone escape probability code of iPoelman fc SpaansI (|2005( ) is used for the radiative transfer in the 
cooling lines (under statistical equilibrium) . Each zone in a PDR was treated this way, and a correction on the cooling 
curve due to line trapping effects can be derived. 

For the warm and dense gas in the simulation clumps, the opacity is concentrated on small physical scales, leading 
to line trapping. For example, cooling lines like [CII] 158 /xm, [01] 63 ^m and CO all have optical depths significantly 
larger than unity across a single PDR interface. Given that this PDR sub-structure is not resolved in the simulations, 
we feel that these opacity effect should be included. 

The dust temperature is computed as part of the PDR model and is often larger than 30 K. The far-infrared radiation 
field produced by this dust is included in the statistical equilibrium of all cooling agents. Particularly water is quite 
sensitive to dust pumping and is an effective heater in the presence of dust grains warmer than 50 K. 

At temperatures above lO'* K the presence of an ionizing component in the impinging radiation field was found to 
suppress the cooling rate, compared to a coUisionally ionized plasma, by a factor of ~ 3 for T = lO'* — 10^ K. This is 
of little concern in our hydrodynamical simulations. Line transfer in the resonantly scattered Lyman a line benefits 
very strongly from dust absorption and re-radiation in the optically thin far-infrared continuum. 




Fig. Al. 
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